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The stationary points (SPs) of a potential energy landscape play a crucial role in understanding 
many of the physical or chemical properties of a given system. Unless they are found analytically, 
there is, however, no efficient method to obtain all the SPs of a given potential. We introduce 
a novel method, called the numerical polynomial homotopy continuation (NPHC) method, which 
numerically finds all the SPs, and is embarrassingly parallelizable. The method requires the non- 
linearity of the potential to be polynomial-like, which is the case for almost all of the potentials 
arising in physical and chemical systems. We also certify the numerically obtained SPs so that they 
are independent of the numerical tolerance used during the computation. It is then straightforward 
to separate out the local and global minima. As a first application, we take the XY model with 
power-law interaction which is shown to have a polynomial-like non-linearity and apply the method. 



Introduction: A Potential energy landscape (PEL) is 
the hyper-surface of some given potential V{x), with 
X — (xi,X2, ...,xn) being the variables (e.g., position co- 
ordinates, fields etc.). Studying the stationary points 
(SPs), defined by the solutions of the system of N equa- 
tions ^^^^^ = 0,i = l,...,iV, of the PEL is crucial in 
learning many physical or chemical properties of the sys- 
tem described by V{x). The SPs are classified accord- 
ing to the number of negative eigenvalues of the Hessian 
matrix H evaluated at each SP: the SPs with no nega- 
tive eigenvalue are called minima and the SPs with at 
least one negative eigenvalue are called saddles. In sta- 
tistical mechanics, for example, the stationary points of 
the PEL have been shown to be directly related to the 
non-analyticity of thermodynamic quantities (i.e., phase 
transitions) in the respective models [1]; the global min- 
imum of a spin glass model is invaluable for studying its 
equilibrium properties. In theoretical chemistry, study- 
ing the properties of the PEL of supercooled liquids and 
glasses has been a very active area of research [2 E], 
specifically, in the study of Kramer's reaction rate the- 
ory for the thermally activated escape from metastable 
states, and the computation of various physical quanti- 
ties like the diffusion constant using the minima of the 
PEL, etc. [2]. In string theory, the SPs of the PEL of 
various supersymmetric potentials correspond to the so- 
called string vacua. A lot of current activities in the 
string phenomenology areas have been focused on devel- 
oping different methods to find these string vacua 

If all SPs are found analytically for the given V{'^), 
then the problem is settled, obviously. But if the analyt- 
ical solutions are intractable, then one has to rely on al- 
ternative methods. Though finding SPs is of the utmost 
importance in so many areas, there do not exist many 
rigorous methods to find SPs, compared to the number 
of methods for minimizing a potential. 

One such method is the gradient-square minimiza- 
tion method in which one minimizes W = |VF(lzf)p 
using some traditional numerical minimization method 



such as Conjugate Gradient, Simulated Annealing, etc. 
[71 [S]. The minima of W, with further restriction that 
W = 0, are the SPs of V{li'). However, there also ex- 
ist minima of W where > 0, and it has been shown 
that the number of these non-SPs grow as the system size 
increases |9j[l^, thus making the method inefficient. 

The Newton-Raphson method (and its sophisticated 
variants) may also be used to solve the system of N non- 
linear equations ^^gj^^ = 0,i = 1, N. Here, an initial 
random guess is fed and then refined to a given numerical 
precision to obtain a solution of the system |H1 E] . How- 
ever, this method still suffers the major drawbacks of the 
gradient-square minimization method, namely, the possi- 
ble existence of large basins of attraction may lead us to 
repeatedly obtain the same SPs and, more importantly, 
no matter how many times we iterate the respective al- 
gorithms with different initial guesses, we are never sure 
if we have found all SPs. 

If the system of stationary equations has polynomial- 
like non-linearity, then the situation is a bit better: in the 
string theory community, for such a system of equations, 
the symbolic methods based on the Grobner basis tech- 
nique are used to solve the system j3HS| which ensure 
that all the SPs are obtained when the computation fin- 
ishes. Roughly speaking, for a given system of multivari- 
ate polynomial equations (which is known to have only 
isolated solutions), the so-called Buchberger Algorithm 
(BA) or its refined variants can compute a new system of 
equations, called a Grobner basis, in which the first equa- 
tion only consists of one of the variables and the subse- 
quent equations consist of increasing number of variables. 
The solutions of the new system remain the same as the 
original system, but the former is easier to solve. This 
method apparently resolves all of the above mentioned 
problems, however, there are a few other problems: the 
BA is known to have suffered from the exponential space 
complexity, i.e., the memory (Random Access Memory) 
required by the machine blows up exponentially with the 
number of variables, equations, terms in each polynomial. 
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etc.). So even for small sized systems, one may not be 
able to compute a Grobner basis. It is also inefficient 
for the systems with irrational coefficients. Furthermore, 
the BA is highly sequential. 

In this Letter, we present the numerical polynomial 
homotopy continuation (NPHC) method which solves 
all the above problems for systems of equations with 
polynomial-like non-linearity. Numerical continuation 
methods have been around for some time [T^ and the 
polynomial homotopy continuation method has also been 
an active area of research |13l IT¥| (see |15l [TF| for the 
earlier account on the NPHC method for various areas 
in theoretical physics). As with the Grobner basis tech- 
nique, the method extensively uses concepts from com- 
plex algebraic geometry and hence we allow the variables 
to take values from the complex space, even if the physi- 
cally important solutions are only the real ones. So long 
as a system of polynomial equations is known to have 
only isolated solutions, the NPHC guarantees that we 
obtain all complex solutions of the system, from which 
all the real solutions can subsequently be filtered out. 
Numerical Polynomial Homotopy Continuation 
Method: Here we introduce the NPHC method to 
solve a system of multivariate polynomial equations. 
Specifically, we consider a system P{x) = which is 
known to have isolated solutions and where P{x) ~ 
( , . . . ■ ^y*-^^ ). Now, the Classical Bezout Theorem 
asserts that for a system of N polynomial equations in N 
variables that is known to have only isolated solutions, 
the maximum number of solutions in is J^^j^ rf^, 
where di is the degree of the ith polynomial. This bound 
is called the classical Bezout bound (CBB). 

Based on the CBB, a homotopy can be constructed 
as H{x,t) = 7(1 - t)Q{x) + tP{x) = 0, where 7 is 
a random complex number. The new system Q{x) — 
(gi(x), . . . , qM{x)), called the start system, is a system of 
polynomial equations with the following properties: (1) 
the solutions of Q{x) = H{x,0) = are known or can 
be easily obtained. The solutions of Q{x) = are called 
start solutions; (2) the number of solutions of Q{x) = 

, 0) = is equal to the CBB of P{x) = 0, (3) the so- 
lution set of (a?, t) = for < < < 1 consists of a finite 
number of smooth paths, each parametrized hy t € [0, 1), 
and (4) every isolated solution of H{x, 1) = P{x) — 
can be reached by some path originating at a solution of 
H{x, 0) = Q{x) = 0. We can then track all the paths 
corresponding to each start solution from t = to t = 1 
and reach (or diverge from) P{x) = = H{x, 1). It is 
rigorously shown that for a generic value of complex 7, all 
paths are regular for t G [0, 1), i.e., there is no singularity 
along the path |13] . By implementing an efficient path 
tracker algorithm, such as the Euler's predictor and New- 
ton's corrector method, all isolated solutions of P{x) = 
can be obtained. We do not delve into the discussion of 
the actual path tracker algorithms used in practice in this 
Letter, except mentioning that in the path tracker algo- 



rithms used in practice, almost all apparent difficulties 
have been resolved, such as tracking singular solutions, 
multiple roots, solutions at infinity, etc.|13|. 

As a trivial illustration, let us take the univariate 
polynomial, P{x) = x^ — 5 = 0, from [TJj. To find its 
solutions, we may for example choose Q{x) = (x"^ — 1) as 
our start system as it satisfies the twin criteria that it has 
the same number of solutions as the CBB of P{x) and is 
also easily solved to obtain the start solutions a; = ±1. 
(Note that, more generally, Q{x) = {x^^ - 1, ... , - 1) 
is the simplest choice as a start system for the multi- 
variate case.) The problem of getting all solutions of 
P{x) = now reduces to simply tracking the solutions of 
H{x, t) — from i = to t = 1 so that the paths begin- 
ning at a; = ±1 lead us to the actual solutions x = ±-\/5. 
Note that, in general, if there are more start solutions 
than actual solutions, then the remaining paths diverge 
as t approaches 1. 

There are several sophisticated computational pack- 
ages such as PHCpack llTj, Bertini [T5] and HOM4PS2 
[T^ which can be used to solve systems of univariate 
and, more importantly, multivariate polynomial equa- 
tions, and are available as freeware. For the sake of com- 
pleteness, the solutions of the above mentioned univariate 
equation are x = ±(2.236067977 + i 10"^°), i.e., ±v^ up 
to the numerical precision. 

Since each path can be tracked independently of all 
others, the NPHC is known as being embarrassingly par- 
allelizable, a feature that makes it very efficient. 

For the multivariate solution is a set of nu- 

merical values of the variables which satisfies each of the 
equations with a given tolerance, Ag^^ 10"^" in our 
set up). Since the variables are allowed to take complex 
values, all the solutions come with real and imaginary 
parts. A solution is a real solution if the imaginary part 
of each of the variables is less than or equal to a given 
tolerance, Ak (^ 10^'' is a robust tolerance for the equa- 
tions we will be dealing with in the next section, below 
which the number of real solutions does not change). All 
these solutions can be further refined with an arbitrary 
precision up to the machine precision. 

The obvious question at this stage would be if the 
number of real solutions depend on Ar. To resolve this 
issue, we use a very recently developed algorithm called 
alphaCertified which is based on the so-called a-theory to 
certify the real non-singular solutions of polynomial sys- 
tems using both exact rational arithmetic and arbitrary 
precision fioating point arithmetic [20 . This is a remark- 
able step, because using alphaCertified we can prove that 
a solution classified as a real solution is actually a real 
solution independent of Ar, and hence these solutions 
are as good as the exact solutions. 

A First Application: As a first application, 
we choose the XY model with long range power- 
law (algebraically decaying) interaction, defined as 
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for our purposes, the normaUzation constant K = 

(W-l) 

(2EjJi ^)"^ and a e [0,oo) [51]. a = re- 
produces the mean-field XY model and a — > cxd re- 
produces the neaxest-neighbour coupling XY model for 
which all the SPs are analytically obtained recently 
[m Ea Eg. We choose a = 0.75 for which the coef- 
ficients take values from R, unlike for example a — 1 
for which the integers are rational numbers. Hence, we 
are already in the domain of problems where the Grob- 
ner basis technique is inefficient. We impose periodic 
boundary condition, i.e., 9k+N = Ok- Spin glass mod- 
els with power-law interaction have gained a huge in- 
terest recently [21]. The chosen model is one of the 
simplest models of this kind with a continuous symme- 
try. The model is also known as the Kuramoto model 
with power-law interaction. The stationary equations 

are ^ = Ejlf ( ^"^'""'"^^'^t""^'""'"-^'^ ) = 0, for 
k^l,!'..,N. To get rid of the global rotation symmetry 
Ok Sk + (f> where g M, we fix one of the angles, say 
6n = 0, and remove the iVth equation out of the system, 
as done in [TSl [23]. Kastner in [21 used a specific class 
of known solutions, called the stationary wave solutions, 
i.e., e'hi^ = where m, n € {!,..., TV}. Below we 

find that there are many more solutions for this model. 

The above system of equations is not apparently a 
system of polynomial equations. But we can transform 
it into one |15l 116] by first using the trigonometric identi- 
ties, sin(0/j — 9k+j) = sin dk cos 9k+j — sin 9k+j cos 9k etc.; 
abbreviating each sin^^ — Sk and cos Ok = c^, in all the 
iV — 1 equations; and finally adding A^ — 1 additional con- 
straint equations as s^. -I- c| — 1 = 0, for A; = 1, . . . , — 1, 
we get a system of 2(A^ — 1) polynomial equations consist- 
ing of 2(A^ — 1) algebraic variables c^s and s^s as below: 

(w-i) 

V~-^ {SkCk+j ^ Sk+jCk + SkCk-j — Sk-jCk) _ Q 
J = l ■' 

sl+cl-l = Q, (1) 

for k — 1, . . . , A — 1. The removal of the global rotation 
symmetry ensures that the system will have only isolated 
solutions. We can now use the NPHC to solve the Eqs. 
([ij for all the s^s and c^s; filter the real solutions out 
using the above mentioned tolerances; and finally get the 
original 6'-variables back by 9k = tan~^(|^) S (— 7r,7r], 
for all /c = 1, . . . , A^ — 1. The results are as follows. 

Firstly, the number of SPs for N = 3,5,7,9,11,13 
are 6,20,168,972,4774,24830, respectively. The left 
panel in Figure [T] shows the number of SPs as a function 
of A^. The SPs for different models with similar sizes 
have been studied using the gradient minimization and 
Newton-Raphson's methods, though the final number of 
SPs is always open to debate [HI Uni , whereas using the 
NPHC method we can find all the SPs with confidence. 



NoofSPs NoDf SPswilliinile^I 




Figure 1: The left panel is the plot for N vs no. of SPs: the 
circles, squares and diamonds represent data points for the 
algebraic decaying XY model, the one-dimensional nearest- 
neighbour XY model for periodic |15l I23| and anti-periodic 
conditions |22| . respectively. The lines are drawn for a guide 
to the eyes. The right panel is the plot for index I vs no of 
SPs with index I. From top to bottom, iV = 13, 11, 9, 7, 5, 3. 
The normalization on / is A'^ — 1 since we have taken 9n ~ 0. 
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Figure 2: Plot for V/N vs index density, starting from top 
left corner A'' = 7, 9, 11, 13, respectively. 

In the gradient minimization method, it is difficult to 
obtain the SPs with higher indices H]. In the NPHC, 
however, all the SPs are treated equally irrespective of 
their indices. The right panel in Figure [T] shows the index 
// (A^ — 1) vs number of SPs for a given A^. Apparently, 
the number of minima grows linearly in this model, which 
is contrary to the conventional wisdom jHlII^. However, 
the number of SPs indeed grows exponentially with in- 
creasing system size as expected jHl • We also observe 
that the global minimum is always the configuration with 
all 6*^ = 0, i = 1, . . . , A — 1, and that at the global min- 
imum, the energy density V/N is zero. In Figure [2] we 
plot V/N vs / for A = 7,9, 11, 13. The apparent non- 
linear relation between these two quantities in the plots 
show a different behaviour to the linear relationship ob- 
served for the Lennard- Jones models. 

Although it is quite difficult to claim any result in 
the thermodynamic limit using such small size systems, 
we note that the minimum value of the |DetH|^/^ is al- 
ways either at the SPs for which all 9i £ {0,7r} for all 
i = 1, . . . , A — 1, or for the SPs for which 9i+i — 9i = 

mod 27r for alH = 1, . . . , A^, with k being an integer. 
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Figure 3: Plot of V/N vs |Det "Hl^/^ starting from upper left 
corner for N = 7,9, 11, 13, respectively. Though these small 
sized systems can in no way be considered as representatives 
of the model in the thermodynamic limit, we still see the 
expected singularity of |Det H]^/^ at V/N 0.5. 

Using a recently spelled out condition for a class of spin- 
glass models to have phase transition, that |Det "Hj^/^ 
evaluated at a class of stationary points tend to go to 
zero at the critical energy in the TV — oo limit, the phase 
transition in our model would occur at this specific class 
of solutions |25]. For = 7, 9, 11, 13, the corresponding 
plots are drawn in Figure [3] where the plot has started 
being filled up by densely populated points around the 
critical energy density V/N ~ 0.5. Hence, we have re- 
produced the results that Kastner worked out from other 
approach in [21] by just studying these small size sys- 
tems. These particular SPs are usually neither minima 
nor maxima. 

It would be interesting to get all the SPs for a € 
(1,2], find out the relevant SPs for the phase transition 
and extrapolate the analysis in the thermodynamic limit 
as suggested in [21]. Verifying a recent conjecture that 
the critical energy density of a class of spin glass models 
can be computed using the Ising-like SPs only (see|5B] for 
the definition of these SPs and the conjecture), will also 
be another important application of the NPHC method. 

In this Letter, we have described and applied a 
novel method, called the numerical homotopy continu- 
ation (NPHC) method, that finds all SPs of a given po- 
tential, provided that the potential has polynomial-like 
non-linearity. As shown in this Letter, even if the given 
potential is not apparently in the polynomial form, its 
stationary equations could be transformed into a polyno- 
mial form by adding suitable constraint equations. So the 
method is very widely applicable. The ability of finding 
all SPs and it being completely parallelizable makes the 
method quite promising to study many areas of theoreti- 
cal physics and chemistry, for example, finding all the SPs 
and hence all the minima of the Lennard- Jones potential 
and its numerous variants; to obtain the string vacua 
of the models for which the symbolic algebraic geome- 



try methods fail due to their algorithmic complexities; to 
study phase transitions in various spin glass models with 
the above mentioned criterion on the hessian determi- 
nant etc. We anticipate that this work will give a thrust 
to the current research in the related areas. 

DM was supported by the U.S. Department of En- 
ergy grant under contract no. DE-FG02-85ER40237 and 
Science Foundation Ireland grant 08/RFP/PHY1462. 
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